# State Minimum Wage and Mental Health Among Children and Adolescents
# Analyses using the National Survey of Children's Health
# N.M. Kavanagh, M. McConnell, N. Slopen
# September 20, 2024

# Please direct questions about this script file to nolankavanagh@fas.harvard.edu.

# Clear R environment
rm(list = ls())

# Load packages
library(here)         # Working directory
library(readstata13)  # Dataset tools
library(tidyverse)    # Analysis tools
library(dplyr)        # Analysis tools
library(psych)        # Analysis tools
library(lfe)          # Modeling tools
library(survey)       # Survey tools
library(ggplot2)      # Graphing tools
library(usmap)        # Graphing tools
library(cowplot)      # Graphing tools
library(scales)       # Graphing tools
library(arsenal)      # Table tools
library(modelsummary) # Table tools
library(kableExtra)   # Table tools
library(cdlTools)     # FIPS tools

##############################################################################
# NSCH dataset preparation
##############################################################################

# Read NSCH datasets into R
nsch_2016_top <- read.dta13("NSCH data/nsch_2016_topical.dta")
nsch_2017_top <- read.dta13("NSCH data/nsch_2017_topical.dta")
nsch_2018_top <- read.dta13("NSCH data/nsch_2018_topical.dta")
nsch_2019_top <- read.dta13("NSCH data/nsch_2019_topical.dta")
nsch_2020_top <- read.dta13("NSCH data/nsch_2020_topical.dta")
nsch_2021_top <- read.dta13("NSCH data/nsch_2021_topical.dta")
nsch_2022_top <- read.dta13("NSCH data/nsch_2022e_topical.dta")

# Treat stratum as character
# Necessary for a smooth merge
nsch_2016_top$stratum <- as.character(nsch_2016_top$stratum)

# Bind all survey years together
nsch <- bind_rows(nsch_2016_top, nsch_2017_top, nsch_2018_top, nsch_2019_top,
                  nsch_2020_top, nsch_2021_top, nsch_2022_top)

# Redefine strata
# Recommended by the NSCH for multi-year analyses
nsch <- nsch %>% mutate(
  strata = case_when(
    stratum %in% c("1")       ~ "1",
    stratum %in% c("2", "2A") ~ "2"
  ))

##############################################################################
# Minimum wage dataset preparation
##############################################################################

# Read minimum wage dataset into R
min_wage_df <- read.csv("Supplemental data/State minimum wages, raw.csv")

# Note: The cleaning steps below produce virtually identical nominal wages for 1980-2019 as this source:
# https://raw.githubusercontent.com/Lislejoem/Minimum-Wage-by-State-1968-to-2020/master/Minimum%20Wage%20Data.csv
# Lisle did not have official data for 2020 onward, so we must clean the data ourselves.

# Copy forward missing years
# Note: These are only relevant for supplemental analyses
# We have complete information for the main study years (2000 onward)
min_wage_df$X1982 <- min_wage_df$X1981
min_wage_df$X1983 <- min_wage_df$X1981
min_wage_df$X1984 <- min_wage_df$X1981
min_wage_df$X1985 <- min_wage_df$X1981
min_wage_df$X1986 <- min_wage_df$X1981
min_wage_df$X1987 <- min_wage_df$X1981
min_wage_df$X1989 <- min_wage_df$X1988
min_wage_df$X1990 <- min_wage_df$X1988
min_wage_df$X1993 <- min_wage_df$X1992
min_wage_df$X1995 <- min_wage_df$X1994
min_wage_df$X1999 <- min_wage_df$X1998

# Convert dataset to long format
min_wage_long <- min_wage_df %>%
  select(order(colnames(min_wage_df))) %>%
  select(State.or.other.jurisdiction, X1980:X2023) %>%
  pivot_longer(cols = X1980:X2023, names_to = "year", values_to = "bls_min_wage")

# Convert state names to FIPS
# Necessary for a smooth merge
min_wage_long$fipsst <- cdlTools::fips(min_wage_long$State.or.other.jurisdiction, to="FIPS")

# Treat year as numeric
# Remove X before each value
min_wage_long$year <- as.numeric(gsub("X", "", min_wage_long$year))

# Generate new column for cleaned minimum wages
min_wage_long$state_min_wage <- min_wage_long$bls_min_wage

# Correct typo (double period) in select rows
min_wage_long$state_min_wage <- gsub("\\.\\.", "\\.", min_wage_long$state_min_wage)

# For ranges, select both values
min_wage_long <- min_wage_long %>% mutate(
  
  # Get values before hyphen or other symbols
  state_min_wage_1 = case_when(
    grepl("[-–/&]", state_min_wage) ~ gsub("[-–/&].*", "", min_wage_long$state_min_wage),
    TRUE ~ gsub("[-–/&].*", "", min_wage_long$state_min_wage)),
  
  # Get values after hyphen or other symbols
  state_min_wage_2 = case_when(
    grepl("[-–/&]", state_min_wage) ~ gsub(".*[-–/&]", "", min_wage_long$state_min_wage),
    TRUE ~ gsub(".*[-–/&]", "", min_wage_long$state_min_wage))
)

# Remove any remaining symbols or letters
# In practice, discard text starting with symbol
min_wage_long$state_min_wage_1 <- gsub("[-–&\\[\\(].*$", "", min_wage_long$state_min_wage_1)
min_wage_long$state_min_wage_2 <- gsub("[-–&\\[\\(].*$", "", min_wage_long$state_min_wage_2)

# Remove leading dollar signs and treat as numeric
min_wage_long$state_min_wage_1 <- as.numeric(gsub("\\$", "", min_wage_long$state_min_wage_1))
min_wage_long$state_min_wage_2 <- as.numeric(gsub("\\$", "", min_wage_long$state_min_wage_2))

# Get upper and lower ends of ranges
min_wage_long$state_min_wage_u <- pmax(min_wage_long$state_min_wage_1, min_wage_long$state_min_wage_2, na.rm = TRUE)
min_wage_long$state_min_wage_l <- pmin(min_wage_long$state_min_wage_1, min_wage_long$state_min_wage_2, na.rm = TRUE)

# Get federal minimum wages
min_wage_fed <- subset(min_wage_long, State.or.other.jurisdiction == "Federal (FLSA)")
min_wage_fed$federal_min_wage <- min_wage_fed$state_min_wage_1
min_wage_fed <- min_wage_fed %>% select(year, federal_min_wage)

# Merge state and federal minimum wages
min_wage_long <- left_join(min_wage_long, min_wage_fed, by="year")

# Get effective minimum wages
# Higher of state or federal minimums
min_wage_long$effective_min_wage_u <- pmax(min_wage_long$state_min_wage_u, min_wage_long$federal_min_wage, na.rm = TRUE)
min_wage_long$effective_min_wage_l <- pmax(min_wage_long$state_min_wage_l, min_wage_long$federal_min_wage, na.rm = TRUE)

# Read CPI dataset into R
# Source: https://www.usinflationcalculator.com/inflation/consumer-price-index-and-annual-percent-changes-from-1913-to-2008/
cpi_df <- read.csv("Supplemental data/CPI by month, raw.csv")

# Select January CPI estimates
# Note: Wage data are effective as of January 1 in each year
cpi_df$year        <- cpi_df$Year
cpi_df$cpi_january <- cpi_df$Jan
cpi_df <- cpi_df %>% select(cpi_january, year)

# Merge CPI estimates
min_wage_long <- left_join(min_wage_long, cpi_df, by="year")

# Compute inflation-adjusted wages
# Use CPI in January 2020: 257.971
min_wage_long$inflation_min_wage_u <- 257.971 / min_wage_long$cpi_january * min_wage_long$effective_min_wage_u
min_wage_long$inflation_min_wage_l <- 257.971 / min_wage_long$cpi_january * min_wage_long$effective_min_wage_l

# Reorder dataset
min_wage_long <- min_wage_long %>% arrange(fipsst, year)

# Generate lagged minimum wages
min_wage_long <- min_wage_long %>%
  group_by(fipsst) %>%
  mutate(
    
    # Not inflation adjusted
    lag_by_1  = lag(effective_min_wage_l, n=1, default=NA),
    lag_by_2  = lag(effective_min_wage_l, n=2, default=NA),
    lag_by_3  = lag(effective_min_wage_l, n=3,  default=NA),
    lag_by_4  = lag(effective_min_wage_l, n=4,  default=NA),
    lag_by_5  = lag(effective_min_wage_l, n=5,  default=NA),
    lag_by_6  = lag(effective_min_wage_l, n=6,  default=NA),
    lag_by_7  = lag(effective_min_wage_l, n=7,  default=NA),
    lag_by_8  = lag(effective_min_wage_l, n=8,  default=NA),
    lag_by_9  = lag(effective_min_wage_l, n=9,  default=NA),
    lag_by_10 = lag(effective_min_wage_l, n=10, default=NA),
    lag_by_11 = lag(effective_min_wage_l, n=11, default=NA),
    lag_by_12 = lag(effective_min_wage_l, n=12, default=NA),
    lag_by_13 = lag(effective_min_wage_l, n=13, default=NA),
    lag_by_14 = lag(effective_min_wage_l, n=14, default=NA),
    lag_by_15 = lag(effective_min_wage_l, n=15, default=NA),
    lag_by_16 = lag(effective_min_wage_l, n=16, default=NA),
    lag_by_17 = lag(effective_min_wage_l, n=17, default=NA),
    
    # Inflation adjusted
    lag_by_1_inf  = lag(inflation_min_wage_l, n=1,  default=NA),
    lag_by_2_inf  = lag(inflation_min_wage_l, n=2,  default=NA),
    lag_by_3_inf  = lag(inflation_min_wage_l, n=3,  default=NA),
    lag_by_4_inf  = lag(inflation_min_wage_l, n=4,  default=NA),
    lag_by_5_inf  = lag(inflation_min_wage_l, n=5,  default=NA),
    lag_by_6_inf  = lag(inflation_min_wage_l, n=6,  default=NA),
    lag_by_7_inf  = lag(inflation_min_wage_l, n=7,  default=NA),
    lag_by_8_inf  = lag(inflation_min_wage_l, n=8,  default=NA),
    lag_by_9_inf  = lag(inflation_min_wage_l, n=9,  default=NA),
    lag_by_10_inf = lag(inflation_min_wage_l, n=10, default=NA),
    lag_by_11_inf = lag(inflation_min_wage_l, n=11, default=NA),
    lag_by_12_inf = lag(inflation_min_wage_l, n=12, default=NA),
    lag_by_13_inf = lag(inflation_min_wage_l, n=13, default=NA),
    lag_by_14_inf = lag(inflation_min_wage_l, n=14, default=NA),
    lag_by_15_inf = lag(inflation_min_wage_l, n=15, default=NA),
    lag_by_16_inf = lag(inflation_min_wage_l, n=16, default=NA),
    lag_by_17_inf = lag(inflation_min_wage_l, n=17, default=NA)
  )

# Merge minimum wage and NSCH data
nsch_all <- left_join(nsch, min_wage_long, by=c("fipsst", "year"))

# Clear old datasets from R
rm(nsch_2016_top, nsch_2017_top, nsch_2018_top, nsch_2019_top,
   nsch_2020_top, nsch_2021_top, nsch_2022_top)

##############################################################################
# Supplemental Medicaid dataset preparation
##############################################################################

# Read Medicaid datasets into R
medicaid_1_5_df  <- read.csv("Supplemental data/Medicaid eligibility, ages 1-5, state-year, raw.csv")
medicaid_6_18_df <- read.csv("Supplemental data/Medicaid eligibility, ages 6-18, state-year, raw.csv")

# Add columns for 2001 and 2007
# Duplicate data from 2000 and 2006
medicaid_1_5_df$X2001  <- medicaid_1_5_df$X2000
medicaid_1_5_df$X2007  <- medicaid_1_5_df$X2006
medicaid_6_18_df$X2001 <- medicaid_6_18_df$X2000
medicaid_6_18_df$X2007 <- medicaid_6_18_df$X2006

# Restructure data to long format
library(reshape2)
medicaid_1_5_df  <- melt(medicaid_1_5_df, id.vars = c("X"),
                         variable.name = "year", value.name = "elig_1_5")
medicaid_6_18_df <- melt(medicaid_6_18_df, id.vars = c("X"),
                         variable.name = "year", value.name = "elig_6_18")

# Clean year variable
# Drop leading "X" in strings
medicaid_1_5_df$year  <- substring(medicaid_1_5_df$year,  2)
medicaid_6_18_df$year <- substring(medicaid_6_18_df$year, 2)

# Rename columns
colnames(medicaid_1_5_df)  <- c("state", "year", "elig_1_5")
colnames(medicaid_6_18_df) <- c("state", "year", "elig_6_18")

# Convert names to FIPS
# Necessary for smooth merge
medicaid_1_5_df$fipsst  <- cdlTools::fips(medicaid_1_5_df$state,  to="FIPS")
medicaid_6_18_df$fipsst <- cdlTools::fips(medicaid_6_18_df$state, to="FIPS")

# Treat year as numeric
medicaid_1_5_df$year  <- as.numeric(medicaid_1_5_df$year)
medicaid_6_18_df$year <- as.numeric(medicaid_6_18_df$year)

# Merge Medicaid and YRBS data
nsch_all <- left_join(nsch_all, medicaid_1_5_df,  by=c("fipsst", "year"))
nsch_all <- left_join(nsch_all, medicaid_6_18_df, by=c("fipsst", "year"))

##############################################################################
# Supplemental EITC dataset preparation
##############################################################################

# Read EITC dataset into R
eitc <- read.csv("Supplemental data/EITC policies, state-year, clean.csv")

# Merge EITC and YRBS data
nsch_all <- left_join(nsch_all, eitc, by=c("fipsst", "year"))

##############################################################################
# Supplemental TANF dataset preparation
##############################################################################

# Read TANF dataset into R
tanf <- read.csv("Supplemental data/TANF benefits family 3, state-year, clean.csv")

# Merge TANF and YRBS data
nsch_all <- left_join(nsch_all, tanf, by=c("fipsst", "year"))

##############################################################################
# Variable preparation
##############################################################################

# Treat fixed effects as factors
nsch_all$year   <- as.factor(nsch_all$year)
nsch_all$fipsst <- as.factor(nsch_all$fipsst)

# Child's age
nsch_all <- nsch_all %>% mutate(
  age = case_when(
    sc_age_years %in% c(0:17) ~ sc_age_years
  ))

# Treat age as factor
nsch_all$age <- as.factor(nsch_all$age)

# Birth year
# Note: For 2019-2022, there is a reported birth year variable
# However, the data quality is questionable for many respondents
# As such, we use the survey year and child's reported age
nsch_all <- nsch_all %>% mutate(
  birth_year = case_when(
    # Subtract age from survey year
    TRUE ~ as.numeric(as.character(year)) - sc_age_years
  ))

# Dichotomize age for interaction models
nsch_all <- nsch_all %>% mutate(
  age_cat = case_when(
    sc_age_years %in% c(13:17) ~ 0, # Adolescents
    sc_age_years %in% c(0:12)  ~ 1  # All other children
  ))

# Child's sex
nsch_all <- nsch_all %>% mutate(
  sex = case_when(
    sc_sex == 1 ~ "Male",
    sc_sex == 2 ~ "Female"
  ))
nsch_all$sex <- factor(nsch_all$sex, levels = c("Male", "Female"))

# Child's race/ethnicity
nsch_all <- nsch_all %>% mutate(
  race_eth = case_when(
    sc_hispanic_r == 1    ~ "Hispanic/Latino",
    sc_race_r == 1        ~ "White",
    sc_race_r == 2        ~ "Black or African American",
    sc_race_r == 3        ~ "American Indian or Alaska Native",
    sc_race_r %in% c(4:5) ~ "Asian, Native Hawaiian, or Pacific Islander",
    sc_race_r %in% c(6:7) ~ "Other or two or more races",
  ))
nsch_all$race_eth <- factor(nsch_all$race_eth, levels = c("American Indian or Alaska Native", "Asian, Native Hawaiian, or Pacific Islander", "Black or African American", "Hispanic/Latino", "White", "Other or two or more races"))

# Dichotomoize race/ethnicity
# Black/Latino vs. other for interaction models
nsch_all <- nsch_all %>% mutate(
  race_eth_cat = case_when(
    sc_hispanic_r == 1      ~ 0, # Black or Hispanic/Latino
    sc_race_r == 2          ~ 0, # Black or Hispanic/Latino
    sc_race_r %in% c(1,3:7) ~ 1  # All other races
  ))

# Adults' highest educational attainment
nsch_all <- nsch_all %>% mutate(
  adult_edu = case_when(
    higrade_tvis == 1 ~ "Less than high school",
    higrade_tvis == 2 ~ "High school (including vocational)",
    higrade_tvis == 3 ~ "Some college or associate degree",
    higrade_tvis == 4 ~ "College degree or higher",
    TRUE ~ "Not provided"
  ))
nsch_all$adult_edu <- factor(nsch_all$adult_edu, levels = c("Less than high school", "High school (including vocational)", "Some college or associate degree", "College degree or higher", "Not provided"))

# Dichotomize educational attainment
# High school (or less) vs. some college (or more)
nsch_all <- nsch_all %>% mutate(
  adult_edu_cat = case_when(
    higrade_tvis %in% c(1:2) ~ 0, # High school or less
    higrade_tvis %in% c(3:4) ~ 1  # All other education levels
  ))

# Generate mean estimated FPL
# Later NSCH years generated 6 imputed FPLs if a household was missing income
nsch_all$fpl_mean <- rowMeans(cbind(nsch_all$fpl_i1, nsch_all$fpl_i2, nsch_all$fpl_i3,
                                    nsch_all$fpl_i4, nsch_all$fpl_i5, nsch_all$fpl_i6), na.rm=T)

# Household's federal poverty level
nsch_all <- nsch_all %>% mutate(
  fpl_category = case_when(
    fpl %in% c(50:99)   | fpl_mean < 100 ~ "Less than 100%",
    fpl %in% c(100:199) | fpl_mean < 200 ~ "100% to 199%",
    fpl %in% c(200:299) | fpl_mean < 300 ~ "200% to 299%",
    fpl %in% c(300:399) | fpl_mean < 400 ~ "300% to 399%",
    fpl %in% c(400:999) | fpl_mean < 999 ~ "400% or greater"
  ))
nsch_all$fpl_category <- factor(nsch_all$fpl_category, levels = c("Less than 100%", "100% to 199%", "200% to 299%", "300% to 399%", "400% or greater"))

# Dichotomize FPL
# Low-income (<200% FPL) vs. higher-income
nsch_all <- nsch_all %>% mutate(
  low_income = case_when(
    fpl %in% c(50:199)  | fpl_mean < 200 ~ 0, # Lower income
    fpl %in% c(200:999) | fpl_mean < 999 ~ 1  # Higher income
  ))

# Family structure
nsch_all <- nsch_all %>% mutate(
  family_struc = case_when(
    family %in% c(1,3)  | family_r %in% c(1,3)  ~ "Two parents, married",
    family %in% c(2,4)  | family_r %in% c(2,4)  ~ "Two parents, not married",
    family %in% c(5:8)  | family_r %in% c(5:6)  ~ "Single parent",
    family %in% c(1:99) | family_r %in% c(1:99) ~ "Another family structure",
    TRUE ~ "Not provided"
  ))
nsch_all$family_struc <- factor(nsch_all$family_struc, levels = c("Two parents, married", "Two parents, not married", "Single parent", "Another family structure", "Not provided"))

# Household nativity
nsch_all <- nsch_all %>% mutate(
  nativity = case_when(
    house_gen == 1 ~ "First-generation household",
    house_gen == 2 ~ "Second-generation household",
    house_gen == 3 ~ "Third-generation household or higher",
    TRUE ~ "Not provided"
  ))
nsch_all$nativity <- factor(nsch_all$nativity, levels = c("First-generation household", "Second-generation household", "Third-generation household or higher", "Not provided"))

# Dichotomize household nativity
# First/second-generation vs. higher
nsch_all <- nsch_all %>% mutate(
  nativity_cat = case_when(
    house_gen %in% c(1:2) ~ 0,
    house_gen %in% c(3)   ~ 1
  ))

# Current depression
nsch_all <- nsch_all %>% mutate(
  depression = case_when(
    k2q32b %in% c(1)                      ~ 1,
    k2q32b %in% c(2) | k2q32a %in% c(1:2) ~ 0
  ))

# Current moderate or severe depression
nsch_all <- nsch_all %>% mutate(
  dep_mod_sev = case_when(
    k2q32c %in% c(2:3)                                         ~ 1,
    k2q32c %in% c(1) | k2q32b %in% c(1:2) | k2q32a %in% c(1:2) ~ 0
  ))

# Current diagnosed anxiety
nsch_all <- nsch_all %>% mutate(
  anxiety = case_when(
    k2q33b %in% c(1)                      ~ 1,
    k2q33b %in% c(2) | k2q33a %in% c(1:2) ~ 0
  ))

# Current moderate or severe anxiety
nsch_all <- nsch_all %>% mutate(
  anx_mod_sev = case_when(
    k2q33c %in% c(2:3)                                         ~ 1,
    k2q33c %in% c(1) | k2q33b %in% c(1:2) | k2q33a %in% c(1:2) ~ 0
  ))

# Current ADD/ADHD
nsch_all <- nsch_all %>% mutate(
  adhd = case_when(
    k2q31b %in% c(1)                      ~ 1,
    k2q31b %in% c(2) | k2q31a %in% c(1:2) ~ 0
  ))

# Current moderate or severe ADD/ADHD
nsch_all <- nsch_all %>% mutate(
  adhd_mod_sev = case_when(
    k2q31c %in% c(2:3)                                         ~ 1,
    k2q31c %in% c(1) | k2q31b %in% c(1:2) | k2q31a %in% c(1:2) ~ 0
  ))

# Current behavior problems
nsch_all <- nsch_all %>% mutate(
  behavior = case_when(
    k2q34b %in% c(1)                      ~ 1,
    k2q34b %in% c(2) | k2q34a %in% c(1:2) ~ 0
  ))

# Current moderate or severe behavior problems
nsch_all <- nsch_all %>% mutate(
  beh_mod_sev = case_when(
    k2q34c %in% c(2:3)                                         ~ 1,
    k2q34c %in% c(1) | k2q34b %in% c(1:2) | k2q34a %in% c(1:2) ~ 0
  ))

# Stomach/digestive issues
nsch_all <- nsch_all %>% mutate(
  stomach_r = case_when(
    stomach %in% c(1) ~ 1,
    stomach %in% c(2) ~ 0
  ))

# Unmet health care (any)
nsch_all <- nsch_all %>% mutate(
  unmet_needs = case_when(
    k4q27 %in% c(1) ~ 1,
    k4q27 %in% c(2) ~ 0,
  ))

# Unmet mental health care
nsch_all <- nsch_all %>% mutate(
  unmet_mental = case_when(
    k4q28x04 %in% c(1)   ~ 1,
    k4q27    %in% c(1:2) ~ 0
  ))

# Missed days of school
# Dichotomize is as 0-6 vs 7+ days
nsch_all <- nsch_all %>% mutate(
  missed_school = case_when(
    k7q02r_r %in% c(4:5)   ~ 1,
    k7q02r_r %in% c(1:3,6) ~ 0
  ))

# Child's employment
nsch_all <- nsch_all %>% mutate(
  child_job = case_when(
    k7q38 == 1 ~ 1,
    k7q38 == 2 ~ 0
  ))

# Lifetime minimum wage
# With and without inflation adjustments
nsch_all <- nsch_all %>% mutate(
  wage_life_nom = case_when(
    sc_age_years == 0 ~ effective_min_wage_l,
    sc_age_years == 1 ~ (effective_min_wage_l + lag_by_1)/2,
    sc_age_years == 2 ~ (effective_min_wage_l + lag_by_1 + lag_by_2)/3,
    sc_age_years == 3 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3)/4,
    sc_age_years == 4 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4)/5,
    sc_age_years == 5 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4 + lag_by_5)/6,
    sc_age_years == 6 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4 + lag_by_5 + lag_by_6)/7,
    sc_age_years == 7 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4 + lag_by_5 + lag_by_6 + lag_by_7)/8,
    sc_age_years == 8 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4 + lag_by_5 + lag_by_6 + lag_by_7 +
                           lag_by_8)/9,
    sc_age_years == 9 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                           lag_by_4 + lag_by_5 + lag_by_6 + lag_by_7 +
                           lag_by_8 + lag_by_9)/10,
    sc_age_years == 10 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4 + lag_by_5 + lag_by_6 + lag_by_7 +
                            lag_by_8 + lag_by_9 + lag_by_10)/11,
    sc_age_years == 11 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4 + lag_by_5 + lag_by_6  + lag_by_7 +
                            lag_by_8 + lag_by_9 + lag_by_10 + lag_by_11)/12,
    sc_age_years == 12 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4 + lag_by_5 + lag_by_6  + lag_by_7 +
                            lag_by_8 + lag_by_9 + lag_by_10 + lag_by_11 +
                            lag_by_12)/13,
    sc_age_years == 13 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4  + lag_by_5 + lag_by_6  + lag_by_7 +
                            lag_by_8  + lag_by_9 + lag_by_10 + lag_by_11 +
                            lag_by_12 + lag_by_13)/14,
    sc_age_years == 14 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 + 
                            lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                            lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                            lag_by_12 + lag_by_13 + lag_by_14)/15,
    sc_age_years == 15 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                            lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                            lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15)/16,
    sc_age_years == 16 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4  + lag_by_5 +  lag_by_6 +  lag_by_7 +
                            lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                            lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15 +
                            lag_by_16)/17,
    sc_age_years == 17 ~ (effective_min_wage_l + lag_by_1 + lag_by_2 + lag_by_3 +
                            lag_by_4  + lag_by_5  + lag_by_6  + lag_by_7 +
                            lag_by_8  + lag_by_9  + lag_by_10 + lag_by_11 +
                            lag_by_12 + lag_by_13 + lag_by_14 + lag_by_15 +
                            lag_by_16 + lag_by_17)/18
  ))
nsch_all <- nsch_all %>% mutate(
  wage_life_inf = case_when(
    sc_age_years == 0 ~ inflation_min_wage_l,
    sc_age_years == 1 ~ (inflation_min_wage_l + lag_by_1_inf)/2,
    sc_age_years == 2 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf)/3,
    sc_age_years == 3 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf)/4,
    sc_age_years == 4 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf)/5,
    sc_age_years == 5 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf + lag_by_5_inf)/6,
    sc_age_years == 6 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf + lag_by_5_inf + lag_by_6_inf)/7,
    sc_age_years == 7 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf + lag_by_5_inf + lag_by_6_inf + lag_by_7_inf)/8,
    sc_age_years == 8 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf + lag_by_5_inf + lag_by_6_inf + lag_by_7_inf +
                           lag_by_8_inf)/9,
    sc_age_years == 9 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                           lag_by_4_inf + lag_by_5_inf + lag_by_6_inf + lag_by_7_inf +
                           lag_by_8_inf + lag_by_9_inf)/10,
    sc_age_years == 10 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf + lag_by_5_inf + lag_by_6_inf + lag_by_7_inf +
                            lag_by_8_inf + lag_by_9_inf + lag_by_10_inf)/11,
    sc_age_years == 11 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf + lag_by_5_inf + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf + lag_by_9_inf + lag_by_10_inf + lag_by_11_inf)/12,
    sc_age_years == 12 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf + lag_by_5_inf + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf + lag_by_9_inf + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf)/13,
    sc_age_years == 13 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf  + lag_by_5_inf + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf  + lag_by_9_inf + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf + lag_by_13_inf)/14,
    sc_age_years == 14 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf + 
                            lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf + lag_by_13_inf + lag_by_14_inf)/15,
    sc_age_years == 15 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf)/16,
    sc_age_years == 16 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf  + lag_by_5_inf +  lag_by_6_inf +  lag_by_7_inf +
                            lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf +
                            lag_by_16_inf)/17,
    sc_age_years == 17 ~ (inflation_min_wage_l + lag_by_1_inf + lag_by_2_inf + lag_by_3_inf +
                            lag_by_4_inf  + lag_by_5_inf  + lag_by_6_inf  + lag_by_7_inf +
                            lag_by_8_inf  + lag_by_9_inf  + lag_by_10_inf + lag_by_11_inf +
                            lag_by_12_inf + lag_by_13_inf + lag_by_14_inf + lag_by_15_inf +
                            lag_by_16_inf + lag_by_17_inf)/18
  ))

# State has EITC program
nsch_all <- nsch_all %>% mutate(
  has_eitc = case_when(
    federal_pct > 0 ~ 1,
    TRUE ~ 0
  ))

# Generate nested sampling clusters
nsch_all$cluster <- paste0(nsch_all$strata, "-", nsch_all$fipsst)

##############################################################################
# Table: Demographic characteristics
##############################################################################

# Code inclusion/exclusion criteria
# That is, must have age and at least once outcome
nsch_all <- nsch_all %>% mutate(
  included_in_study = case_when(
    !is.na(age) & (!is.na(depression) | !is.na(anxiety) | !is.na(adhd) |
      !is.na(behavior) | !is.na(stomach_r) | !is.na(unmet_needs) |
      !is.na(unmet_mental)) ~ 1,
    TRUE ~ 0
  ))

# Assess missingness due to inclusion/exclusion criteria
prop.table(table(nsch_all$included_in_study, useNA="always"))

# Complete cases for models
nsch_all_model <- nsch_all %>%
  
  # Exclude children less than 3
  filter(age %in% c(3:17)) %>%
  
  # Exclude children missing age or at least one outcome
  # Note: All other covariates have no missingness due to "Not provided" options
  filter_at(vars(effective_min_wage_l, included_in_study), all_vars(!is.na(.)))

# Rescale weight so mean is 1
nsch_all_model$weights <- nsch_all_model$fwc/mean(nsch_all_model$fwc, na.rm=T)

# Child characteristics: unweighted
summary(tableby(~ as.numeric(as.character(age)) + sex + race_eth + family_struc + adult_edu + nativity,
                nsch_all_model, digits.pct=1), text=T)

# Child characteristics: weighted
summary(tableby(~ as.numeric(as.character(age)) + sex + race_eth + family_struc + adult_edu + nativity,
                nsch_all_model, digits.pct=1, weights=weights), text=T)

# Mental health outcomes: weighted
summary(tableby(~ factor(depression), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(anxiety), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(adhd), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(behavior), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(stomach_r), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(unmet_needs), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(unmet_mental), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(missed_school), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)
summary(tableby(~ factor(child_job), nsch_all_model, digits.pct=1,
                weights=weights, na.action=na.tableby(TRUE)), text=T)

# Get missingness by outcome
summary(tableby(~ is.na(depression) + is.na(anxiety) + is.na(adhd) + is.na(behavior) +
                  is.na(stomach_r) + is.na(unmet_needs) + is.na(unmet_mental),
                nsch_all_model, digits.pct=1), text=T)
summary(tableby(~ is.na(missed_school) + is.na(child_job),
                subset(nsch_all_model, age %in% c(6:17)), digits.pct=1), text=T)

# Get range of observed minimum wage changes
wage_descriptives <- min_wage_long %>%
  filter(year %in% c(2016:2022), fipsst %in% unique(nsch_all_model$fipsst)) %>%
  group_by(fipsst) %>%
  summarise(min_wage = min(effective_min_wage_l, na.rm=T),
            max_wage = max(effective_min_wage_l, na.rm=T),
            change   = max_wage - min_wage)

# Range of minimum wages
min(wage_descriptives$min_wage); max(wage_descriptives$max_wage)

# Number of states that never changed
table(wage_descriptives$change > 0)

# Range of changes among states that changed
summary(subset(wage_descriptives, change > 0)$change)

##############################################################################
# Graph: Correlations between outcomes
##############################################################################

# Select desired variables
nsch_corr <- nsch_all_model %>% select(depression, anxiety, adhd, behavior, stomach_r,
                                       unmet_needs, unmet_mental, missed_school, child_job)

# Rename columns
nsch_corr <- nsch_corr %>% rename(
  `Depression`         = depression,
  `Anxiety`            = anxiety,
  `ADD/ADHD`           = adhd,
  `Behavioral prob.`   = behavior,
  `Digestive issues`   = stomach_r,
  `Any unmet care`     = unmet_needs,
  `Unmet mental care`  = unmet_mental,
  `7+ school absences` = missed_school,
  `Employment`         = child_job
)

# Run correlation matrix
library(Hmisc); library(corrplot)
corr_matrix <- rcorr(as.matrix(nsch_corr), type="pearson")

# Create matrices for coefficients, p-values
corr_coef <- corr_matrix$r
corr_pval <- corr_matrix$P

# Define colors for correlation plot
col <- colorRampPalette(c("#4477AA", "#77AADD", "#FFFFFF", "#EE9988", "#BB4444"))

# Plot correlation matrix
corr_plot <- corrplot(corr_coef, method = "color", col = col(200),
                      cl.pos = "n", type = "upper", 
                      
                      # Reorder variables
                      order = "original",
                      
                      # Add coefficient of correlation
                      addCoef.col = "black", number.digits = 2,
                      
                      # Labels color, size, and rotation
                      tl.col = "black", tl.cex = 0.7, number.cex = 0.7, tl.srt = 45,
                      
                      # Hide correlation coefficient on principal diagonal
                      diag = FALSE)
print(corr_plot)

##############################################################################
# Functions for TWFE models
##############################################################################

# Function to extract values from TWFE models.
# Requires: Df for saving coefficients, "lfe" model, title of model.
# Returns: Df of values, ready to pass to "clean_coef_df".
make_coef_df <- function(coef_df, model, TITLE) {
  
  # Get outcome and category labels
  if (model$lhs == "depression") {
    outcome  <- "Currently has depression"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "dep_mod_sev") {
    outcome  <- "Currently has moderate or severe depression"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "anxiety") {
    outcome  <- "Currently has anxiety"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "anx_mod_sev") {
    outcome  <- "Currently has moderate or severe anxiety"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "adhd") {
    outcome  <- "Currently has ADD/ADHD"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "adhd_mod_sev") {
    outcome  <- "Currently has moderate or severe ADD/ADHD"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "behavior") {
    outcome  <- "Currently has behavior problems"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "beh_mod_sev") {
    outcome  <- "Currently has moderate or severe behavior problems"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "stomach_r") {
    outcome  <- "Chronic digestive issues in past year"
    category <- "Diagnoses & symptoms"}
  
  if (model$lhs == "unmet_needs") {
    outcome  <- "Any unmet health care in past year"
    category <- "Health care"}
  
  if (model$lhs == "unmet_mental") {
    outcome  <- "Unmet mental health care in past year"
    category <- "Health care"}
  
  if (model$lhs == "missed_school") {
    outcome  <- "7+ school absences in past year"
    category <- "School & work"}
  
  if (model$lhs == "child_job") {
    outcome  <- "Paid employment in past year"
    category <- "School & work"}
  
  # Add row to coefficient df
  coef_df <- rbind(coef_df, cbind(
    outcome, category, TITLE, model$coef[1], model$cse[1], length(model$residuals)))
  
  # Return coefficient df
  return(coef_df)
}

# Function to clean dataframe of TWFE coefficients.
# Requires: Df with columns specified by "make_coef_df".
# Returns: Cleaned dataframe, ready to use in "print_coef_plot".
clean_coef_df <- function(coef_df) {
  
  # Treat as dataframe
  coef_df <- as.data.frame(coef_df)
  
  # Name columns
  colnames(coef_df) <- c("Outcome", "Category", "Sample", "effect", "se", "n")
  
  # Treat columns as numeric
  coef_df$effect <- as.numeric(coef_df$effect)
  coef_df$se     <- as.numeric(coef_df$se)
  coef_df$n      <- as.numeric(coef_df$n)
  
  # Reorder outcomes
  coef_df$Outcome <- factor(
    coef_df$Outcome, levels=c("Paid employment in past year", "7+ school absences in past year", "Unmet mental health care in past year", "Any unmet health care in past year", "Chronic digestive issues in past year", "Currently has moderate or severe behavior problems", "Currently has behavior problems", "Currently has moderate or severe ADD/ADHD", "Currently has ADD/ADHD", "Currently has moderate or severe anxiety", "Currently has anxiety", "Currently has moderate or severe depression", "Currently has depression"))
  
  # Reorder categories
  coef_df$Category <- factor(
    coef_df$Category, levels=c("Diagnoses & symptoms", "Health care", "School & work"))
  
  # Reorder samples
  coef_df$Sample <- factor(
    coef_df$Sample, levels=rev(c(
      
      # Standard models
      "All children (fully adjusted)",
      "All children (FE only)",
      
      # Subgroup models
      "Less than 200% FPL",
      "Adults with high school or less",
      "Black or Hispanic/Latino",
      "First- or second-generation",
      "Adolescents, ages 13–17",
      "Not urban or principal city",
      
      # Robustness checks
      "All children (2020 dollars)",
      "All children (lagged wage)",
      
      # Lifetime wage models
      "All children (fully adj.), nominal wages",
      "All children (FE only), nominal wages",
      "All children (fully adj.), real wages",
      "All children (FE only), real wages",
      
      # Alternate cluster models
      "All children (fully adj., state clusters)",
      "All children (fully adj., nested clusters)",
      "All children (FE only, state clusters)",
      "All children (FE only, nested clusters)"
    )))
  
  # Return df
  return(coef_df)
}

# Function to generate TWFE coefficient plots.
# Requires: Df of coefficients y title, and y limits.
# Returns: Formatted coefficient plot.
print_coef_plot <- function(coef_df, Y_TITLE, Y_MIN, Y_MAX) {
  
  # Generate coefficient plot
  plot <- ggplot(coef_df, aes(x=Outcome, y=effect, group=Sample, color=Sample)) +
    
    # Null line
    geom_hline(yintercept=0, color="black", linewidth=0.25) +
    
    # Point estimates with distinct shapes
    geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
    scale_shape_manual(values = 1:nlevels(coef_df$Sample)) +
    
    # Confidence intervals
    geom_errorbar(aes(ymin=effect-1.96*se, ymax=effect+1.96*se),
                  position = position_dodge(width=0.6), width=0, linewidth=0.8) +
    geom_errorbar(aes(ymin=effect-2.94*se, ymax=effect+2.94*se),
                  position = position_dodge(width=0.6), width=0.3, alpha=0.5) +
    
    # Titles
    ylab(Y_TITLE) +
    ggtitle("All children (ages 3–17), 2016–2022") +
    
    # Theme modifications
    theme_test() +
    theme(legend.position = "bottom",
          text = element_text(size = 10, face = "bold"),
          plot.title = element_text(hjust=0.5),
          axis.title.y = element_blank(),
          axis.ticks = element_blank(),
          strip.background = element_blank(),
          legend.title = element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
          panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
    scale_y_continuous(limits = c(Y_MIN, Y_MAX),
                       breaks = seq(-0.1, 0.1, 0.02),
                       minor_breaks = seq(-0.1, 0.1, 0.01),
                       labels = function(x) paste0(x*100," pp")) +
    scale_color_grey(start=0.7, end=0) +
    facet_grid(Category~., scales="free", space="free_y") +
    coord_flip() +
    guides(shape = guide_legend(reverse=T),
           color = guide_legend(reverse=T))
  
  # Return plot
  return(plot)
}

##############################################################################
# Standard TWFE models
##############################################################################

# Depression
twfe_min_dep_1 <- felm(depression ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_dep_2 <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Anxiety
twfe_min_anx_1 <- felm(anxiety ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_anx_2 <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# ADD/ADHD
twfe_min_add_1 <- felm(adhd ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_add_2 <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Behavior problems
twfe_min_beh_1 <- felm(behavior ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_beh_2 <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Digestive issues
twfe_min_dig_1 <- felm(stomach_r ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_dig_2 <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Any unmet care
twfe_min_unm_1 <- felm(unmet_needs ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_unm_2 <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Unmet mental care
twfe_min_men_1 <- felm(unmet_mental ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_men_2 <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# 7+ school absences
twfe_min_sch_1 <- felm(missed_school ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_sch_2 <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Employment
twfe_min_job_1 <- felm(child_job ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_job_2 <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Get values from models
twfe_df <- NULL

# All children (FE only)
twfe_df <- make_coef_df(twfe_df, twfe_min_dep_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_anx_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_add_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_beh_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_dig_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_unm_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_men_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_sch_1, "All children (FE only)")
twfe_df <- make_coef_df(twfe_df, twfe_min_job_1, "All children (FE only)")

# All children (fully adjusted)
twfe_df <- make_coef_df(twfe_df, twfe_min_dep_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_anx_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_add_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_beh_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_dig_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_unm_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_men_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_sch_2, "All children (fully adjusted)")
twfe_df <- make_coef_df(twfe_df, twfe_min_job_2, "All children (fully adjusted)")

# Clean dataframe of coefficients
twfe_df <- clean_coef_df(twfe_df)

# Get min. and max. N
min(twfe_df$n); max(twfe_df$n)

# Generate coefficient plot: Standard TWFE (for main text)
plot_twfe_1 <- ggplot(twfe_df, aes(x=Outcome, y=effect, group=Sample, color=Sample)) +
  geom_hline(yintercept=0, color="black", linewidth=0.25) +
  geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
  scale_shape_manual(values = 1:nlevels(twfe_df$Sample)) +
  geom_errorbar(aes(ymin=effect-1.96*se, ymax=effect+1.96*se),
                position = position_dodge(width=0.6), width=0.2) +
  ylab("\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes") +
  ggtitle("All children (ages 3–17), 2016–2022") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        plot.title = element_text(hjust=0.5),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(limits = c(-0.045, 0.045),
                     breaks = seq(-0.1, 0.1, 0.02),
                     minor_breaks = seq(-0.1, 0.1, 0.01),
                     labels = function(x) paste0(x*100," pp")) +
  scale_color_grey(start=0.7, end=0) +
  facet_grid(Category~., scales="free", space="free_y") +
  coord_flip() +
  guides(shape = guide_legend(reverse=T),
         color = guide_legend(reverse=T))

# Export figure
ggsave(plot=plot_twfe_1, file="Exhibits/NSCH coefficient plot, TWFE 1.pdf",
       width=6, height=5.5, units='in', dpi=600)

# Generate coefficient plot: Standard TWFE (for appendix)
plot_twfe_2 <- print_coef_plot(
  twfe_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Export figure
ggsave(plot=plot_twfe_2, file="Exhibits/NSCH coefficient plot, TWFE 2.pdf",
       width=6, height=5.5, units='in', dpi=600)

# Variable names for table
var_names <- 
  c("effective_min_wage_l" = "$1 increase in minimum wage",
    "age4"  = "4 years old",
    "age5"  = "5 years old",
    "age6"  = "6 years old",
    "age7"  = "7 years old",
    "age8"  = "8 years old",
    "age9"  = "9 years old",
    "age10" = "10 years old",
    "age11" = "11 years old",
    "age12" = "12 years old",
    "age13" = "13 years old",
    "age14" = "14 years old",
    "age15" = "15 years old",
    "age16" = "16 years old",
    "age17" = "17 years old",
    "sexFemale" = "Female",
    "race_ethAmerican Indian or Alaska Native" =
      "American Indian or Alaska Native",
    "race_ethAsian, Native Hawaiian, or Pacific Islander" =
      "Asian, Native Hawaiian, or Pacific Islander",
    "race_ethBlack or African American" = "Black or African American",
    "race_ethHispanic/Latino" = "Hispanic/Latino",
    "race_ethWhite" = "White",
    "race_ethOther or two or more races" = "Other or two or more races",
    "family_strucTwo parents, not married" = "Two parents, not married",
    "family_strucSingle parent" = "Single parent",
    "family_strucAnother family structure" = "Another family structure",
    "family_strucNot provided" = "Family structure not provided",
    "adult_eduHigh school (including vocational)" = "High school (including vocational)",
    "adult_eduSome college or associate degree" = "Some college or associate degree",
    "adult_eduCollege degree or higher" = "College degree or higher",
    "adult_eduNot provided" = "Adult educatio not provided",
    "nativitySecond-generation household" = "Second-generation household",
    "nativityThird-generation household or higher" = "Third-generation household or higher",
    "nativityNot provided" = "Nativity not provided",
    "elig_1_5"    = "Medicare income eligibility limit, ages 1-5",
    "elig_6_18"   = "Medicare income eligibility limit, ages 6-18",
    "has_eitc"    = "State has EITC",
    "federal_pct" = "State EITC as percent of federal",
    "refundable"  = "State EITC is refundable",
    "max_bft_3"   = "Maximum TANF benefit for family of 3")

# Compile results into tables
modelsummary(list("FE only"    = twfe_min_dep_1,
                  "Fully adj." = twfe_min_dep_2,
                  "FE only"    = twfe_min_anx_1,
                  "Fully adj." = twfe_min_anx_2,
                  "FE only"    = twfe_min_add_1,
                  "Fully adj." = twfe_min_add_2,
                  "FE only"    = twfe_min_beh_1,
                  "Fully adj." = twfe_min_beh_2,
                  "FE only"    = twfe_min_dig_1,
                  "Fully adj." = twfe_min_dig_2,
                  "FE only"    = twfe_min_unm_1,
                  "Fully adj." = twfe_min_unm_2,
                  "FE only"    = twfe_min_men_1,
                  "Fully adj." = twfe_min_men_2,
                  "FE only"    = twfe_min_sch_1,
                  "Fully adj." = twfe_min_sch_2,
                  "FE only"    = twfe_min_job_1,
                  "Fully adj." = twfe_min_job_2),
             gof_omit    = "Log*|AIC|BIC|F|RMSE|Std",
             
             # Options for viewing abbreviated table in R
             coef_omit   = "^(?!effective_min_wage_l)",
             statistic   = c("[{conf.low} to {conf.high}]", "P={p.value}"),
             fmt         = fmt_statistic("estimate"=4, "conf.low"=4, "conf.high"=4, "p.value"=2),
             
             # Options for exporting full table as .csv
             # coef_rename = var_names,
             # statistic   = c("(({std.error}))", "P={p.value}"),
             # fmt         = fmt_statistic("estimate"=4, "std.error"=4, "p.value"=3),
             # output      = "NSCH table, standard TWFE.csv"
             ) %>%
  add_header_above(c(" "=1, "Depression"=2, "Anxiety"=2, "ADD/ADHD"=2,
                     "Behavioral prob."=2, "Digestive issues"=2, "Any unmet care"=2,
                     "Unmet mental care"=2, "7+ school absences"=2, "Employment"=2))

##############################################################################
# TWFE robustness check: Subgroups
##############################################################################

# Generate subgroups of interest
nsch_all_model_l <- subset(nsch_all_model, fpl %in% c(50:199) | fpl_mean < 200) # <200% FPL
nsch_all_model_a <- subset(nsch_all_model, sc_age_years %in% c(13:17)) # Adolescents (ages 13-17)
nsch_all_model_r <- subset(nsch_all_model, sc_hispanic_r == 1 | sc_race_r == 2) # Black or Hispanic/Latino
nsch_all_model_e <- subset(nsch_all_model, higrade_tvis %in% c(1:2)) # High school or less education
nsch_all_model_n <- subset(nsch_all_model, house_gen %in% c(1:2)) # First or second generation
nsch_all_model_u <- subset(nsch_all_model, metro_yn == 2 | mpc_yn == 2) # Not urban or principal city

# Depression
twfe_min_dep_l <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_dep_a <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_dep_r <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_dep_e <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_dep_n <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_dep_u <- felm(depression ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Anxiety
twfe_min_anx_l <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_anx_a <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_anx_r <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_anx_e <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_anx_n <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_anx_u <- felm(anxiety ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# ADD/ADHD
twfe_min_add_l <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_add_a <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_add_r <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_add_e <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_add_n <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_add_u <- felm(adhd ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Behavior problems
twfe_min_beh_l <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_beh_a <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_beh_r <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_beh_e <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_beh_n <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_beh_u <- felm(behavior ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Digestive issues
twfe_min_dig_l <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_dig_a <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_dig_r <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_dig_e <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_dig_n <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_dig_u <- felm(stomach_r ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Any unmet care
twfe_min_unm_l <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_unm_a <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_unm_r <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_unm_e <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_unm_n <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_unm_u <- felm(unmet_needs ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Unmet mental care
twfe_min_men_l <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_men_a <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_men_r <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_men_e <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_men_n <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_men_u <- felm(unmet_mental ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# 7+ school absences
twfe_min_sch_l <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_sch_a <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_sch_r <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_sch_e <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_sch_n <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_sch_u <- felm(missed_school ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Employment
twfe_min_job_l <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_l,
                       weights = nsch_all_model_l$weights)
twfe_min_job_a <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_a,
                       weights = nsch_all_model_a$weights)
twfe_min_job_r <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_r,
                       weights = nsch_all_model_r$weights)
twfe_min_job_e <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_e,
                       weights = nsch_all_model_e$weights)
twfe_min_job_n <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_n,
                       weights = nsch_all_model_n$weights)
twfe_min_job_u <- felm(child_job ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model_u,
                       weights = nsch_all_model_u$weights)

# Get values from models
sub_df <- NULL

# Less than 200% FPL
sub_df <- make_coef_df(sub_df, twfe_min_dep_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_anx_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_add_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_beh_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_dig_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_unm_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_men_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_sch_l, "Less than 200% FPL")
sub_df <- make_coef_df(sub_df, twfe_min_job_l, "Less than 200% FPL")

# Adolescents, ages 13–17
sub_df <- make_coef_df(sub_df, twfe_min_dep_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_anx_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_add_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_beh_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_dig_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_unm_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_men_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_sch_a, "Adolescents, ages 13–17")
sub_df <- make_coef_df(sub_df, twfe_min_job_a, "Adolescents, ages 13–17")

# Black or Hispanic/Latino
sub_df <- make_coef_df(sub_df, twfe_min_dep_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_anx_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_add_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_beh_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_dig_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_unm_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_men_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_sch_r, "Black or Hispanic/Latino")
sub_df <- make_coef_df(sub_df, twfe_min_job_r, "Black or Hispanic/Latino")

# Adults with high school or less
sub_df <- make_coef_df(sub_df, twfe_min_dep_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_anx_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_add_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_beh_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_dig_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_unm_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_men_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_sch_e, "Adults with high school or less")
sub_df <- make_coef_df(sub_df, twfe_min_job_e, "Adults with high school or less")

# First- or second-generation
sub_df <- make_coef_df(sub_df, twfe_min_dep_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_anx_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_add_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_beh_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_dig_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_unm_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_men_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_sch_n, "First- or second-generation")
sub_df <- make_coef_df(sub_df, twfe_min_job_n, "First- or second-generation")

# Not metropolitan principal city
sub_df <- make_coef_df(sub_df, twfe_min_dep_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_anx_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_add_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_beh_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_dig_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_unm_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_men_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_sch_u, "Not urban or principal city")
sub_df <- make_coef_df(sub_df, twfe_min_job_u, "Not urban or principal city")

# Clean dataframe of coefficients
sub_df <- clean_coef_df(sub_df)

# Get min. and max. N
# Before adding standard TWFE models
min(sub_df$n); max(sub_df$n)

# Add standard TWFE models for comparison
sub_df <- rbind(sub_df, twfe_df %>% filter(Sample == "All children (fully adjusted)"))

# Generate coefficient plot: Subgroups
plot_sub <- print_coef_plot(
  sub_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_sub <- plot_sub +
  guides(shape = guide_legend(reverse=T, nrow=4),
         color = guide_legend(reverse=T, nrow=4))

# Export figure
ggsave(plot=plot_sub, file="Exhibits/NSCH coefficient plot, subgroup.pdf",
       width=6, height=8, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Alternate specifications
##############################################################################

# Depression
twfe_min_dep_x <- felm(depression ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_dep_y <- felm(depression ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Anxiety
twfe_min_anx_x <- felm(anxiety ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_anx_y <- felm(anxiety ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# ADD/ADHD
twfe_min_add_x <- felm(adhd ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_add_y <- felm(adhd ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Behavioral problems
twfe_min_beh_x <- felm(behavior ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_beh_y <- felm(behavior ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Digestive issues
twfe_min_dig_x <- felm(stomach_r ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_dig_y <- felm(stomach_r ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Any unmet care
twfe_min_unm_x <- felm(unmet_needs ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_unm_y <- felm(unmet_needs ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Unmet mental care
twfe_min_men_x <- felm(unmet_mental ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_men_y <- felm(unmet_mental ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# 7+ school absences
twfe_min_sch_x <- felm(missed_school ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_sch_y <- felm(missed_school ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Employment
twfe_min_job_x <- felm(child_job ~ inflation_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_job_y <- felm(child_job ~ lag_by_1 +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Get values from models
rob_df <- NULL

# All children (2020 dollars)
rob_df <- make_coef_df(rob_df, twfe_min_dep_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_anx_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_add_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_beh_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_dig_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_unm_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_men_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_sch_x, "All children (2020 dollars)")
rob_df <- make_coef_df(rob_df, twfe_min_job_x, "All children (2020 dollars)")

# All children (lagged wage)
rob_df <- make_coef_df(rob_df, twfe_min_dep_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_anx_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_add_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_beh_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_dig_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_unm_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_men_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_sch_y, "All children (lagged wage)")
rob_df <- make_coef_df(rob_df, twfe_min_job_y, "All children (lagged wage)")

# Clean dataframe of coefficients
rob_df <- clean_coef_df(rob_df)

# Add standard TWFE models for comparison
rob_df <- rbind(rob_df, twfe_df %>% filter(Sample == "All children (fully adjusted)"))

# Get min. and max. N
min(rob_df$n); max(rob_df$n)

# Generate coefficient plot: Alternate specifications
plot_rob <- print_coef_plot(
  rob_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_rob <- plot_rob +
  guides(shape = guide_legend(reverse=T, nrow=3),
         color = guide_legend(reverse=T, nrow=3))

# Export figure
ggsave(plot=plot_rob, file="Exhibits/NSCH coefficient plot, robustness.pdf",
       width=6, height=5.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Logistic regression
##############################################################################

# Treat fixed effects as factors
nsch_all_model$year       <- as.factor(nsch_all_model$year)
nsch_all_model$fipsst     <- as.factor(nsch_all_model$fipsst)
nsch_all_model$birth_year <- as.factor(nsch_all_model$birth_year)

# Define complex sampling design
# Use state clusters and apply survey weights
design_nsch <- svydesign(id=~fipsst, weights=~fwc, data=nsch_all_model)

# Depression
log_min_dep_1 <- svyglm(depression ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_dep_2 <- svyglm(depression ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Anxiety
log_min_anx_1 <- svyglm(anxiety ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_anx_2 <- svyglm(anxiety ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# ADD/ADHD
log_min_add_1 <- svyglm(adhd ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_add_2 <- svyglm(adhd ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Behavioral problems
log_min_beh_1 <- svyglm(behavior ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_beh_2 <- svyglm(behavior ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Digestive issues
log_min_dig_1 <- svyglm(stomach_r ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_dig_2 <- svyglm(stomach_r ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Any unmet care
log_min_unm_1 <- svyglm(unmet_needs ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_unm_2 <- svyglm(unmet_needs ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Unmet mental care
log_min_men_1 <- svyglm(unmet_mental ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_men_2 <- svyglm(unmet_mental ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# 7+ school absences
log_min_sch_1 <- svyglm(missed_school ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_sch_2 <- svyglm(missed_school ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Employment
log_min_job_1 <- svyglm(child_job ~ effective_min_wage_l +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")
log_min_job_2 <- svyglm(child_job ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 +
                          year + fipsst + birth_year,
                        design = design_nsch, family="quasibinomial")

# Get values from models
log_df <- as.data.frame(rbind(
  # Depression
  cbind("Currently has depression", "Diagnoses & symptoms", "All children (FE only)",
        coef(log_min_dep_1)[2],
        SE(log_min_dep_1)[2],
        length(log_min_dep_1$residuals)),
  
  cbind("Currently has depression", "Diagnoses & symptoms", "All children (fully adjusted)",
        coef(log_min_dep_2)[2],
        SE(log_min_dep_2)[2],
        length(log_min_dep_2$residuals)),
  
  # Anxiety
  cbind("Currently has anxiety", "Diagnoses & symptoms", "All children (FE only)",
        coef(log_min_anx_1)[2],
        SE(log_min_anx_1)[2],
        length(log_min_anx_1$residuals)),
  
  cbind("Currently has anxiety", "Diagnoses & symptoms", "All children (fully adjusted)",
        coef(log_min_anx_2)[2],
        SE(log_min_anx_2)[2],
        length(log_min_anx_2$residuals)),
  
  # ADD/ADHD
  cbind("Currently has ADD/ADHD", "Diagnoses & symptoms", "All children (FE only)",
        coef(log_min_add_1)[2],
        SE(log_min_add_1)[2],
        length(log_min_add_1$residuals)),
  
  cbind("Currently has ADD/ADHD", "Diagnoses & symptoms", "All children (fully adjusted)",
        coef(log_min_add_2)[2],
        SE(log_min_add_2)[2],
        length(log_min_add_2$residuals)),
  
  # Behavior problems
  cbind("Currently has behavior problems", "Diagnoses & symptoms", "All children (FE only)",
        coef(log_min_beh_1)[2],
        SE(log_min_beh_1)[2],
        length(log_min_beh_1$residuals)),
  
  cbind("Currently has behavior problems", "Diagnoses & symptoms", "All children (fully adjusted)",
        coef(log_min_beh_2)[2],
        SE(log_min_beh_2)[2],
        length(log_min_beh_2$residuals)),
  
  # Digestive issues
  cbind("Chronic digestive issues in past year", "Diagnoses & symptoms", "All children (FE only)",
        coef(log_min_dig_1)[2],
        SE(log_min_dig_1)[2],
        length(log_min_dig_1$residuals)),
  
  cbind("Chronic digestive issues in past year", "Diagnoses & symptoms", "All children (fully adjusted)",
        coef(log_min_dig_2)[2],
        SE(log_min_dig_2)[2],
        length(log_min_dig_2$residuals)),
  
  # Any unmet care
  cbind("Any unmet health care in past year", "Health care", "All children (FE only)",
        coef(log_min_unm_1)[2],
        SE(log_min_unm_1)[2],
        length(log_min_unm_1$residuals)),
  
  cbind("Any unmet health care in past year", "Health care", "All children (fully adjusted)",
        coef(log_min_unm_2)[2],
        SE(log_min_unm_2)[2],
        length(log_min_unm_2$residuals)),
  
  # Unmet mental care
  cbind("Unmet mental health care in past year", "Health care", "All children (FE only)",
        coef(log_min_men_1)[2],
        SE(log_min_men_1)[2],
        length(log_min_men_1$residuals)),
  
  cbind("Unmet mental health care in past year", "Health care", "All children (fully adjusted)",
        coef(log_min_men_2)[2],
        SE(log_min_men_2)[2],
        length(log_min_men_2$residuals)),
  
  # 7+ school absences
  cbind("7+ school absences in past year", "School & work", "All children (FE only)",
        coef(log_min_sch_1)[2],
        SE(log_min_sch_1)[2],
        length(log_min_sch_1$residuals)),
  
  cbind("7+ school absences in past year", "School & work", "All children (fully adjusted)",
        coef(log_min_sch_2)[2],
        SE(log_min_sch_2)[2],
        length(log_min_sch_2$residuals)),
  
  # Employment
  cbind("Paid employment in past year", "School & work", "All children (FE only)",
        coef(log_min_job_1)[2],
        SE(log_min_job_1)[2],
        length(log_min_job_1$residuals)),
  
  cbind("Paid employment in past year", "School & work", "All children (fully adjusted)",
        coef(log_min_job_2)[2],
        SE(log_min_job_2)[2],
        length(log_min_job_2$residuals))
))

# Name columns
colnames(log_df) <- c("Outcome", "Category", "Sample", "log_odds", "se", "n")

# Reorder outcomes
log_df$Outcome <- factor(
  log_df$Outcome, levels=c("Paid employment in past year", "7+ school absences in past year", "Unmet mental health care in past year", "Any unmet health care in past year", "Chronic digestive issues in past year", "Currently has moderate or severe behavior problems", "Currently has behavior problems", "Currently has moderate or severe ADD/ADHD", "Currently has ADD/ADHD", "Currently has moderate or severe anxiety", "Currently has anxiety", "Currently has moderate or severe depression", "Currently has depression"))

# Reorder categories
log_df$Category <- factor(
  log_df$Category, levels=c("Diagnoses & symptoms", "Health care", "School & work"))

# Reorder samples
log_df$Sample <- factor(log_df$Sample, levels=rev(c("All children (fully adjusted)",
                                                    "All children (FE only)")))

# Treat columns as numeric
log_df$log_odds <- as.numeric(log_df$log_odds)
log_df$se       <- as.numeric(log_df$se)
log_df$n        <- as.numeric(log_df$n)

# Get min. and max. N
min(log_df$n); max(log_df$n)

# Generate coefficient plot: Logistic
plot_log <- ggplot(log_df, aes(x=Outcome, y=exp(log_odds), group=Sample, color=Sample)) +
  geom_hline(yintercept=0, color="black", linewidth=0.25) +
  geom_point(position = position_dodge(width=0.6), size=1, aes(shape=Sample)) +
  scale_shape_manual(values = 1:nlevels(log_df$Sample)) +
  geom_errorbar(aes(ymin=exp(log_odds - 1.96*se), ymax=exp(log_odds + 1.96*se)),
                position = position_dodge(width=0.6), width=0, linewidth=0.8) +
  geom_errorbar(aes(ymin=exp(log_odds - 2.94*se), ymax=exp(log_odds + 2.94*se)),
                position = position_dodge(width=0.6), width=0.3, alpha=0.5) +
  ylab("\nAssociation of $1 increase in minimum wage\nwith odds of mental health outcomes") +
  ggtitle("All children (ages 3–17), 2016–2022") +
  theme_test() +
  theme(legend.position = "bottom",
        text = element_text(size = 10, face = "bold"),
        plot.title = element_text(hjust=0.5),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        legend.title = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(color="light gray", linewidth=0.5),
        panel.grid.minor.x = element_line(color="light gray", linewidth=0.25)) +
  scale_y_continuous(trans = "log10", limits=c(0.25,4),
                     labels=c("0.25", "0.5", "1\nOdds ratio", "2", "4"),
                     breaks=c(0.25, 0.5, 1, 2, 4),
                     minor_breaks = seq(0.25, 10, 0.25)) +
  scale_color_grey(start=0.7, end=0) +
  facet_grid(Category~., scales="free", space="free_y") +
  coord_flip() +
  guides(shape = guide_legend(reverse=T),
         color = guide_legend(reverse=T))

# Export figure
ggsave(plot=plot_log, file="Exhibits/NSCH coefficient plot, logistic.pdf",
       width=6, height=5.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Level-log models
##############################################################################

# Depression
twfe_min_dep_ln_1 <- felm(depression ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_dep_ln_2 <- felm(depression ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Anxiety
twfe_min_anx_ln_1 <- felm(anxiety ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_anx_ln_2 <- felm(anxiety ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# ADD/ADHD
twfe_min_add_ln_1 <- felm(adhd ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_add_ln_2 <- felm(adhd ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Behavior problems
twfe_min_beh_ln_1 <- felm(behavior ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_beh_ln_2 <- felm(behavior ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Digestive issues
twfe_min_dig_ln_1 <- felm(stomach_r ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_dig_ln_2 <- felm(stomach_r ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Any unmet care
twfe_min_unm_ln_1 <- felm(unmet_needs ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_unm_ln_2 <- felm(unmet_needs ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Unmet mental care
twfe_min_men_ln_1 <- felm(unmet_mental ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_men_ln_2 <- felm(unmet_mental ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# 7+ school absences
twfe_min_sch_ln_1 <- felm(missed_school ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_sch_ln_2 <- felm(missed_school ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Employment
twfe_min_job_ln_1 <- felm(child_job ~ log(effective_min_wage_l) |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)
twfe_min_job_ln_2 <- felm(child_job ~ log(effective_min_wage_l) +
                            age + sex + race_eth + family_struc + adult_edu + nativity +
                            elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                            year + fipsst + birth_year | 0 | fipsst,
                          data    = nsch_all_model,
                          weights = nsch_all_model$weights)

# Get values from models
level_log_df <- NULL

# All children (FE only)
level_log_df <- make_coef_df(level_log_df, twfe_min_dep_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_anx_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_add_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_beh_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_dig_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_unm_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_men_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_sch_ln_1, "All children (FE only)")
level_log_df <- make_coef_df(level_log_df, twfe_min_job_ln_1, "All children (FE only)")

# All children (fully adjusted)
level_log_df <- make_coef_df(level_log_df, twfe_min_dep_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_anx_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_add_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_beh_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_dig_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_unm_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_men_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_sch_ln_2, "All children (fully adjusted)")
level_log_df <- make_coef_df(level_log_df, twfe_min_job_ln_2, "All children (fully adjusted)")

# Clean dataframe of coefficients
level_log_df <- clean_coef_df(level_log_df)

# Rescale coefficients for 10% increase in minimum wage
level_log_df$effect <- level_log_df$effect*10*0.01
level_log_df$se     <- level_log_df$se*10*0.01

# Get min. and max. N
min(level_log_df$n); max(level_log_df$n)

# Generate coefficient plot: Alternate specifications
plot_level_log <- print_coef_plot(
  level_log_df,
  Y_TITLE    = "\nAssociation of 10% increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
)

# Export figure
ggsave(plot=plot_level_log, file="Exhibits/NSCH coefficient plot, level-log.pdf",
       width=6, height=5.5, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Nested clusters
##############################################################################

# Depression
twfe_min_dep_1c <- felm(depression ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_dep_2c <- felm(depression ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Anxiety
twfe_min_anx_1c <- felm(anxiety ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_anx_2c <- felm(anxiety ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# ADD/ADHD
twfe_min_add_1c <- felm(adhd ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_add_2c <- felm(adhd ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Behavioral problems
twfe_min_beh_1c <- felm(behavior ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_beh_2c <- felm(behavior ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Digestive issues
twfe_min_dig_1c <- felm(stomach_r ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_dig_2c <- felm(stomach_r ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Any umet needs
twfe_min_unm_1c <- felm(unmet_needs ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_unm_2c <- felm(unmet_needs ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Unmet mental needs
twfe_min_men_1c <- felm(unmet_mental ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_men_2c <- felm(unmet_mental ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# 7+ school absences
twfe_min_sch_1c <- felm(missed_school ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_sch_2c <- felm(missed_school ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Employment
twfe_min_job_1c <- felm(child_job ~ effective_min_wage_l |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
twfe_min_job_2c <- felm(child_job ~ effective_min_wage_l +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | cluster,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Get values from models
clust_df <- NULL

# All children (FE only, nested clusters)
clust_df <- make_coef_df(clust_df, twfe_min_dep_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_anx_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_add_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_beh_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_dig_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_unm_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_men_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_sch_1c, "All children (FE only, nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_job_1c, "All children (FE only, nested clusters)")

# All children (fully adj., nested clusters)
clust_df <- make_coef_df(clust_df, twfe_min_dep_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_anx_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_add_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_beh_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_dig_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_unm_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_men_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_sch_2c, "All children (fully adj., nested clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_job_2c, "All children (fully adj., nested clusters)")

# Add main models for comparison
# All children (FE only, state clusters)
clust_df <- make_coef_df(clust_df, twfe_min_dep_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_anx_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_add_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_beh_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_dig_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_unm_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_men_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_sch_1, "All children (FE only, state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_job_1, "All children (FE only, state clusters)")

# All children (fully adj., state clusters)
clust_df <- make_coef_df(clust_df, twfe_min_dep_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_anx_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_add_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_beh_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_dig_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_unm_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_men_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_sch_2, "All children (fully adj., state clusters)")
clust_df <- make_coef_df(clust_df, twfe_min_job_2, "All children (fully adj., state clusters)")

# Clean dataframe of coefficients
clust_df <- clean_coef_df(clust_df)

# Get min. and max. N
min(clust_df$n); max(clust_df$n)

# Generate coefficient plot: Nested clusters
plot_clust <- print_coef_plot(
  clust_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_clust <- plot_clust +
  guides(shape = guide_legend(reverse=T, nrow=4),
         color = guide_legend(reverse=T, nrow=4))

# Export figure
ggsave(plot=plot_clust, file="Exhibits/NSCH coefficient plot, nested clusters.pdf",
       width=6, height=6, units='in', dpi=600)

##############################################################################
# TWFE robustness check: Symptom severity
##############################################################################

# Depression, moderate or severe
twfe_min_dms_1 <- felm(dep_mod_sev ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_dms_2 <- felm(dep_mod_sev ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Anxiety, moderate or severe
twfe_min_ams_1 <- felm(anx_mod_sev ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_ams_2 <- felm(anx_mod_sev ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# ADD/ADHD, moderate or severe
twfe_min_adm_1 <- felm(adhd_mod_sev ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_adm_2 <- felm(adhd_mod_sev ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Behavior problems, moderate or severe
twfe_min_bms_1 <- felm(beh_mod_sev ~ effective_min_wage_l |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)
twfe_min_bms_2 <- felm(beh_mod_sev ~ effective_min_wage_l +
                         age + sex + race_eth + family_struc + adult_edu + nativity +
                         elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                         year + fipsst + birth_year | 0 | fipsst,
                       data    = nsch_all_model,
                       weights = nsch_all_model$weights)

# Get values from models
severe_df <- NULL

# All children (FE only)
severe_df <- make_coef_df(severe_df, twfe_min_dms_1, "All children (FE only)")
severe_df <- make_coef_df(severe_df, twfe_min_ams_1, "All children (FE only)")
severe_df <- make_coef_df(severe_df, twfe_min_adm_1, "All children (FE only)")
severe_df <- make_coef_df(severe_df, twfe_min_bms_1, "All children (FE only)")

# All children (fully adjusted)
severe_df <- make_coef_df(severe_df, twfe_min_dms_2, "All children (fully adjusted)")
severe_df <- make_coef_df(severe_df, twfe_min_ams_2, "All children (fully adjusted)")
severe_df <- make_coef_df(severe_df, twfe_min_adm_2, "All children (fully adjusted)")
severe_df <- make_coef_df(severe_df, twfe_min_bms_2, "All children (fully adjusted)")

# Clean dataframe of coefficients
severe_df <- clean_coef_df(severe_df)

# Get min. and max. N
min(severe_df$n); max(severe_df$n)

# Generate coefficient plot: Symptom severity
plot_severe <- print_coef_plot(
  severe_df,
  Y_TITLE    = "\nAssociation of $1 increase in minimum wage\nwith rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
)

# Export figure
ggsave(plot=plot_severe, file="Exhibits/NSCH coefficient plot, symptom severity.pdf",
       width=6, height=3, units='in', dpi=600)

##############################################################################
# Lifetime minimum wage models
##############################################################################

# Depression
life_min_dep_1a <- felm(depression ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_dep_2a <- felm(depression ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_dep_1b <- felm(depression ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_dep_2b <- felm(depression ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Anxiety
life_min_anx_1a <- felm(anxiety ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_anx_2a <- felm(anxiety ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_anx_1b <- felm(anxiety ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_anx_2b <- felm(anxiety ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# ADD/ADHD
life_min_add_1a <- felm(adhd ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_add_2a <- felm(adhd ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_add_1b <- felm(adhd ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_add_2b <- felm(adhd ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Behavioral problems
life_min_beh_1a <- felm(behavior ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_beh_2a <- felm(behavior ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_beh_1b <- felm(behavior ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_beh_2b <- felm(behavior ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Digestive issues
life_min_dig_1a <- felm(stomach_r ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_dig_2a <- felm(stomach_r ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_dig_1b <- felm(stomach_r ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_dig_2b <- felm(stomach_r ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Any unmet care
life_min_unm_1a <- felm(unmet_needs ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_unm_2a <- felm(unmet_needs ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_unm_1b <- felm(unmet_needs ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_unm_2b <- felm(unmet_needs ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Unmet mental care
life_min_men_1a <- felm(unmet_mental ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_men_2a <- felm(unmet_mental ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_men_1b <- felm(unmet_mental ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_men_2b <- felm(unmet_mental ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# 7+ school absences
life_min_sch_1a <- felm(missed_school ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_sch_2a <- felm(missed_school ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_sch_1b <- felm(missed_school ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_sch_2b <- felm(missed_school ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Employment
life_min_job_1a <- felm(child_job ~ wage_life_nom |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_job_2a <- felm(child_job ~ wage_life_nom +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

life_min_job_1b <- felm(child_job ~ wage_life_inf |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)
life_min_job_2b <- felm(child_job ~ wage_life_inf +
                          age + sex + race_eth + family_struc + adult_edu + nativity +
                          elig_1_5 + elig_6_18 + has_eitc + federal_pct + refundable + max_bft_3 |
                          year + fipsst + birth_year | 0 | fipsst,
                        data    = nsch_all_model,
                        weights = nsch_all_model$weights)

# Get values from models
life_df <- NULL

# All children (FE only)
life_df <- make_coef_df(life_df, life_min_dep_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_anx_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_add_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_beh_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_dig_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_unm_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_men_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_sch_1a, "All children (FE only), nominal wages")
life_df <- make_coef_df(life_df, life_min_job_1a, "All children (FE only), nominal wages")

life_df <- make_coef_df(life_df, life_min_dep_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_anx_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_add_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_beh_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_dig_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_unm_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_men_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_sch_1b, "All children (FE only), real wages")
life_df <- make_coef_df(life_df, life_min_job_1b, "All children (FE only), real wages")

# All children (fully adjusted)
life_df <- make_coef_df(life_df, life_min_dep_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_anx_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_add_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_beh_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_dig_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_unm_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_men_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_sch_2a, "All children (fully adj.), nominal wages")
life_df <- make_coef_df(life_df, life_min_job_2a, "All children (fully adj.), nominal wages")

life_df <- make_coef_df(life_df, life_min_dep_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_anx_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_add_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_beh_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_dig_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_unm_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_men_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_sch_2b, "All children (fully adj.), real wages")
life_df <- make_coef_df(life_df, life_min_job_2b, "All children (fully adj.), real wages")

# Clean dataframe of coefficients
life_df <- clean_coef_df(life_df)

# Get min. and max. N
min(life_df$n); max(life_df$n)

# Generate coefficient plot: Lifetime wage
plot_life <- print_coef_plot(
  life_df,
  Y_TITLE    = "\nAssociation of $1 increase in lifetime minimum\nwage with rates of mental health outcomes",
  Y_MIN      = -0.045,
  Y_MAX      =  0.045
  )

# Adjust legend
plot_life <- plot_life +
  guides(shape = guide_legend(reverse=T, nrow=4),
         color = guide_legend(reverse=T, nrow=4))

# Export figure
ggsave(plot=plot_life, file="Exhibits/NSCH coefficient plot, lifetime.pdf",
       width=6, height=6, units='in', dpi=600)
